NPS-MA-95-008 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


SOFTWARE  FOR  THE  STAGGERED  AND 
UNSTAGGERED  TURKEL-ZWAS  SCHEMES  FOR  THE 
SHALLOW  WATER  EQUATIONS  ON  THE  SPHERE 

by 

Francis  X.  Giraldo 
Deny  Neta 

November  1995 


Approved  for  public  release;  distribution  is  unlimited. 

Prepared  for:  Naval  Postgraduate  School 
Monterey,  CA  93943 


OTIC 


19960103  121 


I 


NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 


Rear  Admiral  M.  J.  Evans  Richard  Elster 

Superintendent  Provost 

This  report  was  prepared  in  conjunction  with  research  conducted  for  the  Naval  Postgraduate  School  and 
funded  by  the  Naval  Postgraduate  School. 

Reproduction  of  all  or  part  of  this  report  is  authorized. 

This  rqwrt  was  prepared  by: 

/QC  . 

Francis  X.  Giraldo 
NRC  Research  Associate 


Beny  Neta 

Professor  of  Mathematics 


Reviewed  by: 


Released  by: 


PAUL  J.  ^ARTO 


Dean  of  Research 


REPORT  DOCUMENTATION  PAGE 


rorm  Mpprovea 
0MB  No.  0704-0188 


n A?n.nn  Th>  j  P?-  *  average  |  hour  per  response,  including  the  time  for  reviewing  instruaions.  searching  existing  data  sources. 

i  ^  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 

coupon  ^information,  irKluding  for  r^ucing  this  bidden,  to  Washington  Headquarters  Services.  Direaoratefor  information  Operations  and  Reports.  12T?Jefferson 

Davis  Highway.  Suite  1204.  Arlington.  VA  23202«4302.  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduaion  Project  (0704*0 188).  Washington.  DC  20503. 


1.  AGENCY  USE  ONLY  (Leave  blank)  {2.  REPORT  DATE 


November  29.  1995 


4.  TITLE  AND  SUBTITLE  5.  FUNDING  NUMBERS  ~ 

Software  for  the  Staggered  and  Unstaggered  Turkel-Zwas 
Schemes  for  the  Shallow  Water  Equations  on  the  Sphere 

6.  AUTHOR(S)  - 

Francis  X.  Giraldo  and  Beny  Neta 


3.  REPORT  TYPE  AND  DATES  COVERED 

Technical  Report  1  Oct  1995  -  31  Dec  1995 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  AODRESS(ES} 

Naval  Postgraduate  School 
Monterey,  CA  93943-5000 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

NPS-MA-95-008 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  AODRESS(ES) 


10.  SPONSORING /MONITORING 
AGENCY  REPORT  NUMBER 


Naval  Postgraduate  School 
Monterey,  CA  93943 

11.  SUPPLEMENTARY  NOTES  -  '  —  . 

The  views  expressed  in  this  report  are  those  of  the  authors  and  do  not  reflect  the  official  policy  or 
position  of  the  Department  of  Defense  or  the  United  States  Government. 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


Approved  for  public  release;  distribution  is  unlimited. 


13.  ABSTRACT  (Maximum  200  words) 


A  linear  analysis  of  the  shallow  water  equations  in  spherical  coordinates  for  the  Turkel-Zwas^ 
explicit  large  time-step  scheme  was  presented  by  Neta,  Giraldo  and  Navorf  as  well  as  the 
unstaggered*  Turkel-Zwas  scheme  for  the  solution  of  the  shallow  water  equations  on  the  sphere. 


14.  SUBJEa  TERMS 


shallow  water  equations;  finite  differences;  Turkel-Zwas  scheme;  spherical 
coordinates;  staggering. 


15.  NUMBER  OF  PAGES 

40 


16.  PRICE  CODE 


I  5^CVRITY  classification  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRACT 

UlJttWfflED  WiSEmFIED  ummiED 


NSN  7540-01-280-5500 


Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  Z39-18 
298-102 


1.  INTRODUCTION 


In  this  paper  we  present  the  software  developed  for  the  solution  of  the  shallow  water 
equations  in  spherical  coordinates.  The  unstaggered  (original)  Turkel-Zwas  scheme^  and  the 
staggered^  one  are  both  given. 

The  shallow  water  equations  in  spherical  coordinates  are  given  by 
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Here,  /  is  the  Coriolis  parameter  given  by 


(1) 

(2) 

(3) 


/  =  20  sin  9 


(4) 


where  0  is  the  angular  speed  of  the  rotation  of  the  earth,  h  is  the  height  of  the  homogeneous 
atmosphere,  u  and  v  are  the  zonal  and  meridional  wind  components  respectively,  9  and  A 
are  the  latitudinal  and  longitudinal  directions  respectively,  a  is  the  radius  of  the  earth,  and 
g  is  the  gravitational  constant. 

In  section  2  we  present  the  unstaggered  scheme  (modified  as  suggested  by  Neta^).  In 
section  3  we  present  the  staggered  method  as  developed  by  Neta,  Giraldo  and  Navon^.  In 
section  4  we  present  the  input  file  required  including  a  logical  parameter  to  choose  between 
the  staggered  and  unstaggered  versions.  In  section  5  we  present  the  code  developed. 
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2.  UNSTAGGERED  TURKEL-ZWAS  SCHEME 

The  Turkel-Zwas  scheme  for  the  nonlinear  shallow  water  equations  in  spherical  coordi¬ 
nates  takes  the  following  form: 


+2At  1^(1  -  a){fj  +  t&n$j)vij 
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+f  (/i  +  tan 

Ku  -  ^Lij)  +  vij 

+?  fe+g  -  ^ii-g)] 

^/i+g  H  ^  tan^j+g^ 

+f  (fj-q  +  tan 

-o- +  ^k,j  (^li+i  “  ^fc,i-i) 

[(1  -  «)  (“i+pj  -  4-pj) 


^*2  (“fc+p,i+9  ^fc-p,i+g  "t"  “fc+p,i-g  ^*-p>i-g)]  p 
[(1  -  «)  (^li+g  COS  Oj+g  -  vl-_^  COS  ^,_g) 
(^fc+P,i+g  CO®  ^i+g  ^fc+p,i-g  cos  Oj-q 
+t^fc-pj+g  cos  Oj^q  —  Vk-pJ-q  COS  j  ^  j" 


4 


where 


At  At 
^  a  AX  a  AO' 

For  0  =  1  the  geostrophic  balance  and  the  incompressibility  condition  are  satisfied  to  a 
higher  order  in  the  Cartesian  coordinate  case  (See  Turkel  and  Zwas/  Navon  and  de  Villiers®). 

Note  that  there  is  a  typo  in  equation  (11a)  of  Turkel- Zwas^  which  is  our  equation 
(5).  We  have  also  modified  (to  get  a  symmetric  approximation  as  suggested  by  Neta^  for  a 
rectangular  domain)  the  right  hand  side  of  (11c)  in  Turkel-Zwas^  which  is  (7)  here. 
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3.  STAGGERED  TURKEL-ZWAS  SCHEME 

The  staggered  version  of  the  Turkel-Zwas  scheme  as  proposed  by  Neta,  Giraldo  and 
Navon^  takes  the  following  form: 

4?  =  nt/  -<T 
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where  a  is  given  by  (8). 


4.  INPUT 


The  input  file  contains  four  lines.  The  first  input  line  contains  2  integers: 
nx  =  number  of  longitudinal  points 
ny  =  number  of  latitudinal  points. 

The  second  one  contains  3  integers: 
dt  =  time  step  in  seconds 
time^jjg^]  =  final  time  in  hours 
iplot  =  number  of  iterations  per  plot 

The  third  input  line  contains  2  integers  and  a  real  number: 

( 

p  =  stencil  in  longitudinal  direction 
q  =  stencil  in  latitudinal  direction 
alf=Pade-type  differencing  weighting  factor 

The  last  input  line  contains  2  logical  variables: 
pstag  =  staggering  in  p  if  .true, 
qstag  =  staggering  in  q  if  .true. 

For  example: 

64  32 

100  24  100000 

1  1  0.0 

.false,  .false. 
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5.  CODE 


♦These  lines  of  code  contain  the  pcirameter  statements  for  the  global 
♦definitions  of  many  important  parameters. 

♦ - 

implicit  real^8(a-h,o-z) 
parameter  (  imax=128,  jmax=64  ) 

parameter  (  mx=imax^jmax,  mxpoi=mx,  mxele=mx,  mxbou=mx/5,  nd=4  ) 
parameter  (  tol=1.0e-6,  g=10.0,  rk=0.1  ) 


♦This  program  solves  the  Shallow  Water  Equations 

♦on  a  sphere  with  Periodic  B.C.'s  in  the  latitudinal  direction  (theta) 
♦cind  longitudinal  direction  (lambda)  using  a 
♦Staggered  Turkel-Zwas  Scheme  as  suggested  by  B.  Net  a. 

♦Derivatives  are  obtained  via  2nd  order  differencing  with  some  matching 
♦conditions  developed  by  F.X.  Giraldo  to  satisfy  continuous  derivatives 
♦across  the  poles . 

♦Written  by  F.X.  Giraldo  on  10/95 

♦  NRC  Fellow 

♦  Department  of  Mathematics 

♦  Naval  Postgraduate  School 

♦  Monterey,  CA  93940 


program  nturkel 
include  ’param.h’ 


real  taray(2) 


! global  matrices 
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dimension  f (mxpoi) 
dimension  coord (mxpoi, 2) 
integer  node(imax, jmax) ,  p,  q 
logical  pstag,  qstag 

{♦♦♦primitive  variables  arrays*** 

!u  velocity  arrays 

dimension  urn (mxpoi) ,  uO(mxpoi),  up(mxpoi),  ui (mxpoi) 

!v  velocity  arrays 

dimension  vm (mxpoi) ,  vO (mxpoi),  vp (mxpoi),  vi (mxpoi) 

!phi  arrays 

dimension  phim(mxpoi) ,  phiO(mxpoi),  phip(mxpoi),  phii(mxpoi) 

!Read  the  Input  Variables  and  create  the  Grid 
call  init (phiO ,u0 , vO , phii , ui , vi , node , coord , f , 

$  iipoin,xmin,xmax,ymin,3rmax,comega,nx,ny,dx,dy,dt, 

$  nt ime , rade , iplo t , omega , alpha , velmax , cf 1 , p , q , alf , 

$  pstag, qstag) 

! Calculate  Total  Available  Potential  Energy 
call  energy (ae,aeO,phiO,uO,vO,npoin, time) 
write(*,’("  Energy  =  ",el2.4)')ae 

aeO=ae 

time=0.0 
pi=4. 0*atan(l . 0) 
open ( 1 , f ile= ' phi . out ’ ) 
open(2,file=’u.out ') 
open(3,file='v.out’) 
if  (mod(ntime,iplot) .eq.O)  then 
isets=ntime/iplot  +  1 
else 

isets=ntime/iplot  +  2 
end  if 

write(l,*)isets 


10 


write(2,*)isets 
write(3,*) isets 

call  output (phiO,uO,vO,npoin,tiinejnx,ny,phi_meaii) 

!*+*TIME  MARCH 
t ime l=dt ime (t aray ) 

!Do  itime=l  Eulerian  steps 

do  itime=l jUtime 
time=time  +  dt 
ttime=time/ (3600.0) 

writeC^j'C  timestep  time  (hours)  =  "  ,i5,2x,el2.4)  Oitime,ttime 
if  (itime.eq.l)  then 

call  mat suno (phim , phi 0 , phip ,um,u0,up,vm,v0,vp, 

$  coord , f , npoin , dt , dx , dy , node , nx , ny , r ade , comega , 

$  alpha, p,q,alf) 

else 

call  tzst ag (phim , phiO , phip , um , uO , up , vm , vO , vp , coord , 

$  f , npo in , dt , dx , dy , node , nx , ny , r ade , comega , 

$  alpha, p,q,alf,pstag,qstag) 

endif 

call  sf ilt er (phip , up , vp , node , nx , ny , dx , dy ) 

call  t ime_f ilt er (phim , phi 0 , phip ,  um , uO , up ,  vm , vO , vp , npoin , 

$  it ime) 

call  update(phim,phiO,phip,im,uO,up,vm,vO,vp, npoin) 
if  (mod(itime,iplot) .eq.O) 

$  call  output (phi0,u0,v0, npoin, time, nx,ny,phi_mean) 

call  energy (ae , aeO , phiO , uO , vO , npoin ,t ime) 
write(*,'("  Energy  =  ",el2. 4) ')ae 
end  do 

t ime2=et ime (t aray ) 
tclock= (taray (1) +taray (2) ) 

write(*,'("  Total  CPU  time  in  seconds  =  " ,el2 .4) Otclock 
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! Check  time  for  printing  output 
if  (modCntime.iplot) .ne.O) 

$  call  output (phiO,uO,vO,npo in, time, nx,ny,phi_mean) 
close(l) 

! Compute  the  L2  Error  Norm 

cal 1  norm (phi 0 , uO , vO , phi i , ui , vi , node , coord , dx , dy , nx , ny , 

$  phi_norm,u_norm) 

print*,'  L2  NORM  =  ' ,phi_norm,u_norm 
print*,’  dt  dx  dy  velmax  =  ’ ,dt ,dx,dy,velmax 
print* , ’  **  CFL  =  ’ , cf 1 
stop 
end 

♦ - - 

*This  subroutine  calculates  the  Available  Energy  of  the  2D  Shallow  Water 
*Equations  in  spherical  coordinates 
*Written  by  F.X.  Giraldo  on  10/95 

* - * 

subrout ine  energy ( ae , aeO ,phi,u,v,npoin,t ime ) 
include  ’param.h’ 

! global  arrays 

dimension  phi (mxpoi) ,  u(mxpoi) ,  v(mxpoi) 
ae=0 . 0 

lloop  thru  the  elements 

do  ip=l,npoin 

vel2=u(ip)**2  +  v(ip)**2 
ae=ae  +  (phi(ip))*vel2  +  phi(ip)**2 
end  do 

ae=ae/ (2 . 0*g) 
if  (time.gt .0.0)  then 

if  (ae.gt .l.l*ae0. or. ae.lt. 0.9*ae0)  then 
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write(*,'("  ♦Fatal  Error+  10’/.  Init  Energy  Exceeded! ")  0 

write(+,'("  Current .Energy  Initial_Energy=  " 

$  ,2(lx,el2.4))Oae,ae0 

endif 
endif 
return 
end 

* - * 

♦This  subroutine  reads  in  the  input  file. 

♦The  info  read  is:  the  number  of  grid  points  in  x  eind  y  (nx,ny), 

♦  time  step,  final  time,  and  time  steps  per  plotting, 

♦  P>  alpha 

♦  pstag,  qstag. 

♦where  .true,  means  that  it  is  staggered  and  .false,  means  it  is  unstaggered. 


♦Written  by  F.X.  Giraldo  on  10/95 

♦ - - 

subroutine  init (phiO , uO , vO , phii , ui , vi , node , coord , f , 

$  npo in , xmin , xmax , ymin , ymax , comega , nx , ny , dx , dy , dt , 

$  nt ime , rade , iplot , omega , alpha , velmax , cf 1 , p , q , alf , 

$  pstag, qstag) 


include  'param.h’ 
dimension  coord(mxpoi,2) 

dimension  phiO(mxpoi),  uO(mxpoi),  vO(mxpoi) 
dimension  phii (mxpoi) ,  ui(mxpoi),  vi(mxpoi),  f(mxpoi) 
integer  node(imax, jmax) ,  p,  q 
logical  pstag, qstag 

!Read  Input  File 

read(^,^)nx,ny 
read (♦ , ♦)dt , t ime.f inal , iplot 
read (♦,^)p,q, alf 
read(^,^)pstag, qstag 

! check  bounds 
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if  (nx.gt .imax.or.ny.gt. jmax)  then 

write(*,'("  Error!  -  Need  to  Enlarge  IMAX  and  JMAX")') 
write(*,'("  nx  ny  imax  jmax  =  ",4(i3,lx)) ')nx,ny,imax, jmax 
stop 
end  if 

ISet  some  constants 

pi=4 . 0*atan(l . 0) 
rade=6 . 37e+06 

time_f inal=time_f inal*3600 . 0 

ntime=nint (time.f inal/dt ) 

xmin=0 . 0 

xmax®2 . 0*pi 

ymin=-pi/2.0 

ymax=pi/2.0 

xl=xmax-xmin 

yl=symax-3min 

dx=xl/ (nx) 

dy=yl/(ny) 

phi_mean=5 . 768e4 

omega*20 . 0 

comega=7 . 292e-05 

velmax=-le5 

alpha.f cor=0 . 0 

alpha=0 . 0 

!set  the  Initial  Conditions 

ip=0 

do  j=l,ny 

olat=ymin  +  real(j-0.5)*dy 
do  i=l,nx 

olon=xmin  +  real(i-0.5)*dx 

ip=ip+l 

node(i, j)=ip 
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coord (ip, l)=olon 
coord (ip , 2) =olat 

f (ip)=2.0*comega*(  ~cos(olon)*cos(olat)*sin(alpha_f cor)  + 
$  sin(olat)*cos(alpha_fcor)  ) 

u0(ip)=omega*sin(olon)*(sin(olat)**3  - 
$  3*sin(olat)*cos(olat)**2) 

v0(ip)=omega*sin(olat)**2*cos(olon) 
phiO(ip)=phi_mean  + 

$  2*comega*rade*omega*sin(olat)**3*cos(olat)*sin(olon) 

phii (ip) =phiO (ip) 
ui(ip)=uO(ip) 
vi(ip)=vO(ip) 

vell=abs(uO(ip) )  +  abs(vO(ip))  +  sqrt (2*phi0(ip)) 
velmELX=inax  ( velrnaj: ,  vel  1 ) 
end  do 
end  do 

dl=sqrt (dx**2  +  dy**2) 
cf l=dt*velmax/ (dl*rade) 

print*,'  dt  dx  dy  velmax  =  ' , dt , dx , dy , velmax 

print*,’  **  CFL  =  ’,cfl 

npoin=nx*ny 


return 

end 


*This  subroutine  solves  the  2D  Shallow  Water  Equations  in  Spherical 
*Coordinates  using  a  Staggered  Turkel-Zwas  Scheme. 

*Written  by  F.X.  Giraldo  on  10/95 

* - 

subrout ine  mat  suno (phim , phiO , phip , um , uO , up , vm , vO , vp , 

$  coord, f ,npoin,dt,dx,dy,node,nx,ny,rade,comega, 

$  alpha, p,q,alf) 
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include  ’param.h’ 

dimension  phim(mxpoi) ,  phiO(mxpoi),  pliip(mxpoi) 
dimension  um(mxpoi) ,  uO(mxpoi),  up(mxpoi) 
dimension  vm(mxpoi) ,  vO(mxpoi),  vp(mxpoi) 
dimension  coord(mxpoi,2) ,  f(mxpoi) 
integer  node (imax,j max) ,  p,  q,  ph,  qh 

!Loop  through  the  points  and  integrate  using  Forward  Time 
land  Centered  Space... 

! Predictor  Stage  (forward  Euler) 

ph=p 

qh=q 

nxh=nx/2 

do  i=l,nx  ILoop  through  Longitudinal  Nodes 

il=i-l 
i2=i+l 
i3=i-p 
i4=i+p 
i3h=i-ph 
i4h=i+ph 

! Longitudinal  Periodicity 
if  (il.lt.l)  il=il  +  nx 
if  (i2.gt.nx)  i2=i2  -  nx 

! Longitudinal  Periodicity  -P’s  and  +P’s 
if  (i3.lt. 1)  i3=i3  +  nx 
if  (i4.gt.nx)  i4=i4  -  nx 

! Longitudinal  Periodicity  -P/2’s  and  +P/2’s 
if  (i3h.lt.l)  i3hs=i3h  +  nx 
if  (i4h.gt.nx)  i4h=i4h  -  nx 

ILoop  through  Latitudinal  Nodes 


do  j=l,ny 
j2=j+l 
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j3=j-q 

j4=j+q 

j3h=j-qh 

j4h=j +qh 

jlsign=l 

j2sign=l 

j3sign=l 

j4sign=l 

j3hsign=l 

j4hsign=l 

! South  Pole  Periodicity 

ijl=i 

if  (jl.lt.l)  then 

jl=l 

jlsign=-l 
ijl=ijl  +  nxh 

if  (ijl.gt.nx)  ijl=ijl  -  nx 
endif 

! North  Pole  Periodicity 

ij2=i 

if  Cj2.gt.ny)  then 
j2=ny 
j2sign=-l 
ij2=ij2  +  nxh 

if  (ij2.gt.nx)  ij2=ij2  -  nx 
endif 

! South  Pole  Periodicity  -Q’s 

ij3=i 

ippj3=i4 

impj3=i3 

if  (j3.lt. 1)  then 
j3=l  -  j  +  q 
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j3sign=-l 
ij3=ij3  +  nxh 
ippj3=ippj3  +  nxh 
impj3*inipj3  +  nxh 
if  (ij3-gt.nx)  ij3=ij3  -  nx 
if  (ippj3.gt.nx)  ippj3=ippj3  -  nx 
if  (inipj3.gt .nx)  impj3=impj3  -  nx 
endif 

! North  Pole  Periodicity  +Q's 

ij4=i 

ippj4=i4 

impj4=i3 

if  (j4.gt.ny)  then 
j4=2*ny  +  1  -  j  -  q 
j4sign=-l 
ij4=ij4  +  nxh 
ippj4=ippj4  +  nxh 
impj4=impj4  +  nxh 
if  (ij4.gt.nx)  ij4=ij4  -  nx 
if  (ippj4.gt .nx)  ippj4=ippj4  -  nx 
if  (impj4.gt .nx)  impj4=impj4  -  nx 
endif 

! South  Pole  Periodicity  -Q/2’s 

ij3h=i 
ippj 3h=i4h 
impj3h=i3h 
if  (j3h.lt. 1)  then 
j3h=l  -  j  +  qh 
j3hsign=-l 
ij3h=ij3h  +  nxh 
ippj3h=ippj3h  +  nxh 
impj3h=impj3h  +  nxh 
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if  (ij3h.gt.nx)  ij3h=ij3h  -  nx 
if  (ippj3h.gt .nx)  ippj3h=ippj3h  -  nx 
if  (impj3h.gt .nx)  impj3h=impj3h  -  nx 
endif 

! North  Pole  Periodicity  +Q/2's 

ij4h=i 
ippj4h=i4h 
impj4h=i3h 
if  (j4h.gt.ny)  then 

j4h=2*ny  +  1  -  j  -  qh 
j4hsign=-l 
ij4h=ij4h  +  nxh 
ippj4h=ippj4h  +  nxh 
impj4h*impj4h  +  nxh 
if  (ij4h.gt.nx)  ij4h=ij4h  -  nx 
if  (ippj4h.gt .nx)  ippj4h=ippj4h  -  nx 
if  (impj4h.gt .nx)  impj4h=impj4h  -  nx 
endif 

!Set  up  the  Node  Pointers 

! Centered  Diff  Grid  Points 

ip=node(i, j) 
ipl=node(il,j) 
ip2=node(i2, j) 
jpl=node(ijl,jl) 
jp2=node(ij2,j2) 

!Turkel-Zwas  Grid  Points 

ip3*node(i3, j) 

ip4=node(i4,j) 

jp3=node(ij3,j3) 

jp4=node(ij4,j4) 

ip3 j  p3«node (imp j  3 , j  3) 

ip4 j  p3=node (ipp j  3 , j  3) 
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ip3 j  p4=node ( imp j  4 , j  4) 
ip4 j  p4=node ( ipp j  4 , j  4) 

! Staggered  Grid  Points 

ip3h=node ( i3h , j ) 
ip4h=node (i4h , j ) 
j  p3h=node ( i j  3h , j  3h) 
j p4h=no de ( i j 4h , j 4h ) 
ip3h j  p3h=node ( imp j  3h , j  3h) 
ip4h j  p3h=node ( ipp j  3h , j  3h) 
ip3h j  p4h*node ( imp j  4h , j  4h) 
ip4h j  p4h=node ( ipp j  4h , j  4h) 

! Longitudes  and  Latitudes 

o 1 on=co ord ( ip , 1 ) 
olat=coord (ip , 2) 
olonpp=olon  +  p*dx 
olonmp=olon  -  p*dx 
olonpq=olon 

if  (j4sign.eq.-l)  olonpq=olonpq  +  pi 

olatpq=olat  +  q*dy 

olonmq=olon 

if  (j3sign.eq.-l)  olonmqs'olonmq  +  pi 
olatmq=olat  -  q*dy 

! Staggered  Longitudes  and  Latitudes 

olonpqh=olon 

if  (j4hsign.eq.-l)  olonpqh=olonpqh  +  pi 

olatpqh=olat  +  qh*dy 

olonmqh=olon 

if  (j3hsign.eq.-l)  olonmqh=olonmqh  +  pi 
olatmqh=olat  -  qh*dy 

! Coriolis  Force 

fip=2*comega*(  -cos(olon)*cos(olat)*sin(alpha)  + 

$  sin(olat)*cos (alpha)  ) 
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f ip4=2*comega*(  -cos(olonpp)*cos(olat)*sin(alpha)  + 

$  sin (olat)*cos (alpha)  ) 

f ip3=2*comega*(  ~cos(oloniBp)*cos(olat)*sin(alpha)  + 

$  sin(olat)*cos(alpha)  ) 

f jp4=2*comega*(  ~cos(olonpq)*cos(olatpq)*sin(alpha)  + 

$  sin (olatpq)*cos (alpha)  ) 

f jp3*2*coiiiega*(  -cos (oloniiiq)*cos(olatmq)+sin(alpha)  + 

$  sin (olatmq) *003 (alpha)  ) 

fip=f (ip) 
fip4=f (ip4) 
fip3=f(ip3) 
f  jp4=f  (jp4) 
f jp3=f (jp3) 

! integrate  PHI 
phip(ip)*phiO(ip) 

$  -  dt*uO(ip)/(rade*cos(olat))*(phiO(ip2)-phiO(ipl))/(2*dx) 

$  -  dt*v0(ip)/(rade)*(phi0(jp2)-phi0(jpl))/(2*dy) 

$  -  dt*phiO(ip)/(rade*cos(olat))*( 

^  (1.0-alf)*(  (u0(ip4h)-u0(ip3h))/(2*ph*dz)  + 

$  (j4hsign*vO(jp4h)*cos(olatpqh)  - 

$  j3hsign+v0(jp3h)*cos(olatmqh))/(2*qh*dy)  )  + 

$alf/2*(  (j4hsign*uO(ip4hjp4h)-j4hsign*uO(ip3hjp4h))/(2*ph*dx)  + 

$  (j3hsign*u0(ip4hjp3h)-j3hsign*u0(ip3hjp3h))/(2*ph*dx)  + 

$  (j4hsign*vO(ip4hjp4h)*cos(olatpqh)  - 

$  j3hsign*vO(ip4hjp3h)*cos(olatmqh))/(2*qh*dy)  + 

$  (j4hsign*vO(ip3hjp4h)*cos(olatpqh)  - 

$  j3hsign*vO(ip3hjp3h)*cos(olatmqh))/(2*qh*dy)  )  ) 

phip ( ip ) =phi 0 ( ip ) 

! integrate  U 

up(ip)=uO(ip) 

$  -  dt*uO(ip)/(rade*cos(olat))*(uO(ip2)-uO(ipl))/(2*dx) 
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$  -  dt*v0(ip)/rade=i=(j2sign*u0(jp2)-jlsign*u0(jpl))/(2*dy) 

$  -  dt/(rade*cos(olat))*(phiO(ip4h)-phiO(ip3h))/(2*ph*dx) 

$  +  dt*( 

$  (1.0-alf)*(fip  +  uO(ip)/rade*tan(olat))*vO(ip)  + 

$  alf/2*(fip4  +  uO(ip4)/rade*taii(olat))*vO(ip4)  + 

^  alf/2*(fip3  +  u0(ip3)/rade*tan(olat))*v0(ip3)  ) 


up(ip)=uO(ip) 

! integrate  V 

vp(ip)=vO(ip) 

$  -  dt*u0(ip)/(rade*cos(olat))*(v0(ip2)-v0(ipl))/(2*dx) 

$  -  dt*v0(ip)/rade*(j2sign*v0(jp2)-jlsign*v0(jpl))/(2*dy) 

$  -  dt/rade*(  phiO(jp4h)-phiO(jp3h)  )/(2*qh+dy) 

$  -  dt*C 

$  (l-0-alf)*(fip  +  uO(ip)/rade*taii(olat))*uO(ip)  + 

$  alf/2*(fjp4  + 

^  j^®^6®*'^®(jp4)/rade*tan(olatpq))*j4sign*u0(jp4)  + 

$  alf/2*(fjp3  + 

$  j3sign*u0(jp3)/rade*tan(olatmq))*j3sign*u0(jp3)  ) 


vp(ip)=vO(ip) 
end  do 
end  do 
return 
end 


♦ - 

*This  subroutine  computes  the  L2  Norm 
*for  the  Geopotential  and  Velocity  using 
Trapezoid  Rule  Integration. 

♦Written  by  F.X.  Giraldo  on  10/95 

- - 


subroutine  norm (phi 0, uO ,v0, phi i ,ui 


,  vi , node , coord , dx , dy , nx , ny , 


* 


* 
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$ 


phi.norm , u.norm) 
include  ’param.h' 

dimension  phiO(mxpoi),  uO(mxpoi),  vO(mxpoi) 
dimension  phii(mxpoi) ,  ui(mxpoi),  vi(mxpoi) 
dimension  coord(mxpoi,2) ,  phih(128,64) ,  uh(128,64),  vh(128,64) 
integer  node (imax,j max) 

pi=4 . 0*atan(l . 0) 
open ( 10 ,f ile= ' phih . out ’ ) 
open(20,file=’uh.out ' ) 
open (30 ,f ile= ’ vh . out ' ) 
do  j=l,64 

do  i=l,128 

read(10,*)phih(i, j) 
read (20,*)uh(i,j) 
read(30,*)vh(i, j) 
end  do 
end  do 
close(lO) 
close(20) 
close(30) 
do  j=l,ny 
do  i=l,nx 

ip=node(i, j) 
ui(ip)=uli(2*i-l,2*j-l) 
vi ( ip) »vh ( 2*i- 1 , 2* j -1 ) 
phii(ip)=phih(2=t'i-l,2*j-l) 
end  do 
end  do 

phi_top*0.0 
phi_bot=0 . 0 
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u_top=0.0 
u_bot=0 . 0 
do  j=l,ny-l 
do  i=l,nx-l 
il=node(i, j) 
i2=node(i+l,j) 
i3=node(i+l,j+l) 
i4=node(i, j+1) 
olatl=coord(il,2) 
olat2=coord(i2,2) 
olat3=coord(i3,2) 
olat4=coord(i4,2) 

phil=(phiO(il)  -  phii(il))*cos(olatl) 

ul=(uO(il)  -  ui(il))*cos(olatl) 

vl=(vO(il)  -  vi(il))*cos(olatl) 

phi2=(phi0(i2)  -  phii(i2))*cos(olat2) 

u2=(u0(i2)  -  ui(i2))=i'cos(olat2) 

v2=(v0(i2)  -  vi(i2))*cos(olat2) 

phi3=(phi0(i3)  -  phii(i3))*cos(olat3) 

u3=(u0(i3)  -  ui(i3))*cos(olat3) 

v3=(v0(i3)  -  vi(i3))*cos(olat3) 

phi4=(phi0(i4)  -  phii(i4))*cos(olat4) 

u4=(u0(i4)  -  ui(i4))*cos(olat4) 

v4=(v0(i4)  -  vi(i4))*cos(olat4) 

phi=dx*dy*(phil  +  phi2  +  phi3  +  plii4)/4 

phie=dx*dy*(phii(il)  +  phii(i2)  +  phii(i3)  +  phii(i4))/4 

u=dx*dy*(ul  +  u2  +  u3  +  u4)/4 

ue=dx*dy*(ui(il)  +  ui(i2)  +  ui(i3)  +  ui(i4))/4 

v=dx*dy*(vl  +  v2  +  v3  +  v4)/4 

ve=dx*dy*(vi(il)  +  vi(i2)  +  vi(i3)  +  vi(i4))/4 

phi_top=phi_top  +  (  phi  )**2 
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phi_bot=phi_bot  +  (  phie  )**2 

u_top=u_top  +  (  u  ')**2  +  (  V  )**2 
u_bot=u_bot  +  (  ue  )**2  +  (  ve  )**2 
end  do 

end  do 

phi_nonn=l . 0/ (4 . 0*pi) *sqrt (phi.t op/phi _bot) 
u_norm=l .  0/(4. 0=i'pi)*sqrt  (u_top/u_bot) 

return 

end 

♦This  subroutine  writes  the  output.  It  is  currently  set  only  to 
♦print  the  geopotential  and  wind  velocities  at  each  node  point. 

♦Written  by  F.X.  Giraldo  on  10/95 

* - - 

subrout ine  output (phi , u , v , npoin , t ime , nx , ny , phi.mean) 
include  'param.h' 

dimension  phi(mxpoi),  u(mxpoi),  v(inxpoi) 

pi='4 .  O^ataind .  0) 
dtime=time/3600 . 0 

writed,  >  (2(i6,l2),el6.8)  ')nx,ny,dtime 
writed,'(el2.4)  ')(phi(ip),  ip=l, npoin) 
write(2, ’ (2(i6,lx),el6.8) Onx,ny,dtime 
write(2, ' (el2.4) ') (u(ip) ,  ip=l, npoin) 
write(3, ' (2(i6,lx) ,el6.8) dnx,ny,dtime 
write(3, '(el2.4) ')(v(ip),  ip=l, npoin) 
return 
end 


♦This  subroutine  performs  the  Robert  time  filtering  using  a 
♦Laplacian  type  time-diffusion  term  that  smoothens  the  values  spatially. 
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♦Written  by  F.X.  Giraldo  on  10/95 


subrout ine  sf i It er (phiO , uO , vO , node , nx , ny , dx , dy ) 
include  'param.h' 

dimension  phiO (mxpoi) ,  phip(mxpoi) 
dimension  uO(mxpoi),  up(mxpoi) 
dimension  vO(mxpoi),  vp(mxpoi) 
integer  node (imax,j max) 

do  i=l,nx 
il=i-l 
i2=i+l 

if  (il.lt. 1)  il=nx 
if  (i2.gt.nx)  i2=l 
do  j=l,ny 

if  (j .gt.2.or.j .It .ny-1)  goto  100 
j2»j+l 

!Set  up  the  Node  Pointers 

ip=node(i, j) 
ipl=node(il, j) 
ip2=node(i2, j) 
jpl=node(i,jl) 
jp2=node(i,j2) 

phi0_xx=(  phi0(ip2)  -  2*phi0(ip)  +  phiO(ipl)  )/(dx*dx) 
u0_xx=(  u0(ip2)  -  2*u0(ip)  +  uO(ipl)  )/(dx*dx) 
v0_xx=(  v0(ip2)  -  2*v0(ip)  +  vO(ipl)  )/(dx*dx) 
phi0_yy=(  phi0(jp2)  -  2*phi0(ip)  +  phiO(jpl)  )/(dy*dy) 
u0_yy=(  u0(jp2)  -  2*u0(ip)  +  uO(jpl)  )/(dy*dy) 
v0_yy=(  v0(jp2)  -  2*v0(ip)  +  vO(jpl)  )/(dy*dy) 
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100 


! South  Pole  Periodicity 
if  (jl.lt.l)  then 

jl=l 

ijl=i  +  nx/2 

if  (ijl.gt.nx)  ijl=ijl  -  nx 
jpl=node(ijl,jl) 

phi0_yy=(  phi0(jp2)  -  2*phi0(ip)  +  phiO(jpl)  )/(dy*dy) 
u0_yy=(  u0(jp2)  -  2*u0(ip)  +  uO(jpl)  )/(dy*dy) 
v0_yy=(  v0(jp2)  -  2*v0(ip)  +  vO(jpl)  )/(dy*dy) 
endif 

INorth  Pole  Periodicity 
if  (j2.gt.ny)  then 
j2=ny 

ij2=i  +  nx/2 

if  (ij2,gt.nx)  ij2=ij2  -  nx 
jp2=node(ij2,j2) 

phi0_yy=(  phi0(jp2)  -  2*phi0(ip)  +  phiO(jpl)  )/(dy*dy) 
u0_yy=(  u0(jp2)  -  2*u0(ip)  +  uO(jpl)  )/(dy*dy) 
v0.yy=(  v0(jp2)  -  2*v0(ip)  +  vO(jpl)  )/(dy*dy) 
endif 

phip(ip)=phiO(ip)  +  rk*(  phiO.xx  +  phi0_yy  ) 
up(ip)=u0(ip)  +  rk*(  u0_xx  +  u0_yy  ) 
vp(ip)=v0(ip)  +  rk*(  v0_xx  +  v0_yy  ) 
continue 
end  do 
end  do 

do  i=l,nx 
do  j=l,ny 

if  (j .gt .2.or.j .lt.ny-1)  goto  200 
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ip=node(i, j) 
phiO(ip)=phip(ip) 
uO(ip)=up(ip) 
vO(ip)=vp(ip) 

200  continue 
end  do 
end  do 

return 

end 

* - - 

♦This  subroutine  performs  the  Robert  time  filtering  using  a 

♦Laplacian  type  time-diffusion  term  that  smoothens  the  values  temporally. 

♦Written  by  F.X.  Giraldo  on  10/95 

♦ - - 

subrout ine  t ime_f ilt er (phim , phiO , phip , um , uO , up , vm , vO , vp , npo in , 

$  itime) 

include  ’peLram.h’ 

dimension  phim(mxpoi) ,  phiO(mxpoi),  phip(mxpoi) 
dimension  um(mxpoi) ,  uO(mxpoi),  up(mxpoi) 
dimension  vm(mxpoi) ,  vO(mxpoi),  vp(mxpoi) 

if  (itime. eq. 2)  then 
do  ip=l,npoin 

phi0(ip)=phi0(ip)  +  rk^(  phip(ip)  -  phiO(ip)  ) 
u0(ip)=u0(ip)  +  rk^(  up(ip)  -  uO(ip)  ) 
v0(ip)=v0Cip)  +  rk^(  vp(ip)  -  vO(ip)  ) 
end  do 

else  if  (itime. gt. 2)  then 
do  ip=l,npoin 

phi0(ip)=phi0(ip)  +  rk^(  phip(ip)  -  2^phi0(ip)  +  phim(ip)  ) 
u0(ip)=u0(ip)  +  rk^(  up(ip)  -  2^u0(ip)  +  um(ip)  ) 
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vO(ip)=vO(ip)  +  rk*(  vp(ip)  -  2*v0(ip)  +  vm(ip)  ) 
end  do 
endif 
return 
end 


♦This  subroutine  solves  the  2D  Shallow  Water  Equations  in  Spherical 
♦Coordinates  using  a  Staggered  Turkel-Zwas  Scheme. 

♦Written  by  F.X.  Giraldo  on  10/95 


subrout ine  t zst  ag (phim , phiO , phip ,  im , uO , up , vm , vO , vp , coord , 

$  f , npo in , dt , dx , dy , node , nx , ny , rade , comega , 

$  alpha,p,q,alf ,pstag,qstag) 

include  ’paaram.h’ 

dimension  phim (mxpoi) ,  phiO(mxpoi),  phip(mxpoi) 
dimension  um(mxpoi),  uO(mxpoi),  up(mxpoi) 
dimension  vm (mxpoi) ,  vO(mxpoi),  vp(mxpoi) 
dimension  coord (mxpoi, 2) ,  f (mxpoi) 
integer  node (imaxjjmax) ,  p,  q,  ph,  qh 
logical  pstag,  qstag 

ILoop  through  the  points  and  integrate  using  Forward  Time 
!euid  Centered  Space... 

! Predictor  Stage  (forward  Euler) 

if  (pstag)  then 
ph=p/2 
else 
ph=p 
endif 

if  (qstag)  then 
qh=q/2 
else 
qh=q 
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endif 

alfh=0 

alfu=0 

alfv=0 

nxh=nx/2 

do  i=l,iix  !Loop  through  Longitudinal  Nodes 

il=i-l 
i2=i+l 
i3=i-p 
i4=i+p 
i3h=i-ph 
i4h=i+ph 

! Longitudinal  Periodicity 
if  (il.lt. 1)  il=il  +  nx 
if  (i2.gt.nx)  i2=i2  -  nx 

! Longitudinal  Periodicity  -P’s  and  +P’s 
if  (i3.1t.l)  i3=i3  +  nx 
if  (i4.gt.nx)  i4=i4  -  nx 

! Longitudinal  Periodicity  -P/2’s  and  +P/2's 
if  (i3h.lt. 1)  i3h=i3h  +  nx 
if  (i4h.gt.nx)  i4h=i4h  -  nx 

!Loop  through  Latitudinal  Nodes 

do  j=l,ny 

j2=j+l 

j3=j-q 

j4=j+q 

j3h=j-qh 

j4h=j+qh 

jlsign=l 

j2sign=l 

j3sign=l 
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j4sign=l 

j3hsign=l 

j4hsign=l 

! South  Pole  Periodicity 

ij  l=i 

if  (jl.lt.l)  then 

jl=l 

jlsign=-l 
ijl=ijl  +  nxh 

if  (ijl.gt.nx)  ijl=ijl  -  nx 
endif 

! North  Pole  Periodicity 

ij2=i 

if  (j2.gt.ny)  then 
j2=ny 
j2sign=-l 
ij2=ij2  +  nxh 

if  (ij2.gt.nx)  ij2=ij2  -  nx 
endif 

! South  Pole  Periodicity  -Q's 

ij3=i 
ippj3=i4 
imp j  3=i3 

if  (j3.1t.l)  then 
j3=l  -  j  +  q 
j3sign=-l 
ij3=ij3  +  nxh 
ippj3=ippj3  +  nxh 
impj3=impj3  +  nxh 
if  (ij3.gt.nx)  ij3=ij3  -  nx 
if  (ippj3.gt .nx)  ippj3=ippj3  -  nx 
if  (impj3.gt .nx)  impj3=impj3  -  nx 
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endif 


! North  Pole  Periodicity  +Q’s 

ij4=i 

ippj4=i4 

impj4*i3 

if  (j4.gt.ny)  then 
j4=2*ny  +  1  -  j  -  q 
j4sign=-l 
ij4=ij4  +  nxh 
ippj4=ippj4  +  nxh 
impj4=impj4  +  nxh 
if  (ij4.gt.nx)  ij4=ij4  -  nx 
if  (ippj4.gt .nx)  ippj4=ippj4  -  nx 
if  (impj4.gt.nx)  impj4=impj4  -  nx 
endif 

! South  Pole  Periodicity  -Q/2’s 

ij3h=i 
ippj3h=i4h 
impj 3h=i3h 
if  (j3h.lt. 1)  then 
j3h=l  -  j  +  qh 
j3hsign=-l 
ij3h=ij3h  +  nxh 
ippj3h=ippj3h  +  nxh 
impj 3h= imp j3h  +  nxh 
if  (ij3h.gt.nx)  ij3h=ij3h  -  nx 
if  (ippj3h.gt .nx)  ippj3h=ippj3h  -  nx 
if  (impj3h.gt .nx)  impj3h=impj3h  -  nx 
endif 

! North  Pole  Periodicity  +Q/2's 

ij4h=i 

ippj4h=i4h 
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impj4h=i3h 
if  (j4h.gt.ny)  then 

j4h=2*ny  +  1  -  j  -  qh 
j4hsign=-l 
ij4h=ij4h  +  nxh 
ippj4h=ippj4h  +  nxh 
impj4h=impj4h  +  nxh 
if  (ij4h.gt.nx)  ij4h=ij4h  -  nx 
if  (ippj4h.gt.nx)  ippj4h=ippj4h  -  nx 
if  (impj4h.gt .nx)  impj4h=iiiipj4h  -  nx 
endif 

!Set  up  the  Node  Pointers 

{Centered  Diff  Grid  Points 

ip=node(i, j) 
ipl=node(il,j) 
ip2=node(i2, j) 
jpl=node(ijl,jl) 
jp2=node(ij2,j2) 

!Turkel-Zwas  Grid  Points 

ip3=node(i3,j) 
ip4=node(i4,j) 
jp3=node(ij3,j3) 
jp4=node(ij4, j4) 
ip3 j  p3=node (imp j  3 , j  3) 
ip4 j  p3*node (ipp j  3 , j  3) 
ip3 j  p4=node (imp j  4 , j  4) 
ip4 j  p4=node (ipp j  4 , j  4) 

{Staggered  Grid  Points 

ip3h=node(i3h, j ) 
ip4h=no  de ( i4h , j ) 
jp3h=node(ij3h, j3h) 
jp4h=node(ij4h, j4h) 
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ip3hj  p3h=node (imp j  3h , j  3h) 
ip4h j  p3h=node ( ipp j  3h , j  3h) 
ip3hjp4h=node (impj4h, j4h) 
ip4h j p4h=node ( ipp j 4h , j 4h ) 

! Longitudes  emd  Latitudes 

olon=coord(ip,l) 
olat=coord(ip,2) 
olatpq=olat  +  q*dy 
olatmq=olat  -  q*dy 

! Staggered  Longitudes  and  Latitudes 
olatpqh=olat  +  qh*dy 
olatmqh=olat  -  qh*dy 

! Coriolis  Force 

fip=f (ip) 
fip4=f (ip4) 
f ip3=f (ip3) 
fjp4=f(jp4) 
f jp3=f (jp3) 

! integrate  PHI 
phip ( ip) =phim (ip) 

$  -  dt*u0(ip)/(rade*cos(olat))*(phi0(ip2)-phi0(ipl))/dx 

$  -  dt*v0(ip)/(rade)*(phi0(jp2)-phi0(jpl))/dy 

$  -  dt*phiO(ip)/(rade*cos(olat))*( 

$  (1.0-alf)*(  (u0(ip4h)-u0(ip3h))/(ph*dx)  + 

$  (j4hsign*v0(jp4h)*cos(olatpqh)  - 

$  j3hsign*v0(jp3h)*cos(olatmqh))/(qh*dy)  )  + 

$  alf/2*(  (j4hsign*u0(ip4hjp4h)-j4hsign*u0(ip3hjp4h))/(ph*dx)  + 

$  (j3hsign*u0(ip4hjp3h)-j3hsign*u0(ip3hjp3h))/(ph*dx)  + 

$  (j4hsign*v0(ip4hjp4h)*cos(olatpqh)  - 

$  j3hsign*v0(ip4hjp3h)*cos(olatmqh))/(qh*dy)  + 

$  (j4hsign*v0(ip3hjp4h)*cos(olatpqh)  - 

$  j3hsign*v0(ip3hjp3h)*cos(olatmqh))/(qh+dy)  )  ) 
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c  phip(ip)=phiO(ip) 

! integrate  U 

up(ip)=um(ip) 

$  -  dt*u0(ip)/(rade*cos(olat))*(u0(ip2)-u0(ipl))/dx 

$  -  dt*v0(ip)/rade*(j2sign*u0(jp2)-jlsign*u0(jpl))/dy 

$  -  dt/(rade*cos(olat))*(phiO(ip4h)-phiO(ip3h))/(ph*dx) 

$  +  2*dt*( 

$  (1 .0-alf)*(fip  +  uO(ip)/rade*taai(olat))+vO(ip)  + 

$  alf/2*(fip4  +  uO(ip4)/rade*tan(olat))*vO(ip4)  + 

$  alf/2*(fip3  +  uO(ip3)/rade*tan(olat))*vO(ip3)  ) 

c  up(ip)=uO(ip) 

! integrate  V 

vp(ip)=vm(ip) 

$  -  dt*u0(ip)/(rade+cos(olat))*(v0(ip2)-v0(ipl))/cix 

$  -  dt*v0(ip)/rade*(j2sign*v0(jp2)-jlsign*v0(jpl))/dy 

$  -  dt/rade*(  phiO(jp4h)-phiO(jp3h)  )/(qh*dy) 

$  -  2=t=dt*( 

$  (1.0-alfv)*(fip 

$  +  uO(ip)/rade*taii(olat))*uO(ip) 

$  +  alfv/2*(fjp4 

$  +  j4sign*uO(jp4)/rade*tan(olatpq))*j4sign*uO(jp4) 

$  +  alfv/2*(fjp3 

$  +  j3sign*uO(jp3)/rade*tan(olatmq))*j3sign*uO(jp3)  ) 

c  vp(ip)=vO(ip) 

end  do 
end  do 
return 
end 

* - - 

♦This  subroutine  updates  the  arrays  PHIM,UM,VM,PHI0,U0,V0, 

♦Written  by  F.X.  Giraldo  on  10/95 

♦ - - 


35 


subroutine  updat e ( phim , ph i 0 , ph ip , um , uO , up , vm , vO , vp , npo in) 
include  'param.h' 

dimension  phim (mxpoi) ,  phiO(mxpoi),  phip(mxpoi) 
dimension  um(mxpoi) ,  uO(mxpoi),  up(mxpoi) 
dimension  vm (mxpoi ) ,  vO(mxpoi),  vp(mxpoi) 

!Loop  through  all  the  nodes  and  update 

do  ip=l,npoin 

! Update  F (x-2*alpha, t-dt )=F (x-alpha,t) 

phim (ip) =phiO (ip) 

um(ip)=uO(ip) 

vm(ip)=vO(ip) 

lUpdate  F(x-alpha,t)=F(x,t+dt) 

phiO(ip)=phip(ip) 
uO(ip)=up(ip) 
vO(ip)=vp(ip) 
end  do 

return 

end 
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